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Abstract 

Two-dimensional Ising models on the honeycomb lattice and the square lattice with striped 
random impurities are studied to obtain their phase diagrams. Assuming bimodal distributions 
of the random impurities where all the non-zero couplings have the same magnitude, exact 
critical values for the fraction p of ferromagnetic bonds at the zero-temperature (T = 0) are 
obtained. The critical lines in the p—T plane are drawn by numerically evaluating the Lyapunov 
exponent of random matrix products. 


1 Introduction 

Ising models have attracted a particular attention in the statistical physics as simplest models that 
exhibit phase transitions and critical phenomena. The most popular is the pure Ising ferromagnet 
model on the square lattice, for which not only the transition temperature m but also the free 
energy have been obtained exactly [161 HI]- The hexagonal (honeycomb) lattice is the secondly 
simplest two-dimensional lattice. The pure Ising model on the honeycomb lattice has been studied 
by Wannier [55] and Houtapel m and exactly solved. Apart from the Ising model, quantum spin 
models such as the Kitaev model m or Kitaev-Heisenberg model [ana on the honeycomb lattice 
have recently received a lot of theoretical and experimental attentions [n la [la i5a isni • such a 
growing interest in the honeycomb-lattice systems motivates us to revisit the Ising model on the 
honeycomb lattice. 

In contrast to pure models, less is understood for disordered Ising models. The exact free energy 
has not been obtained in fully random Ising models in two dimension. Although the transition 
temperature and/or phase diagram have been investigated using the replica trick, few is revealed 
so far. Domany showed for the diluted Ising model on the square lattice with concentration p of 
the ferromagnetic bonds that the zero-temperature phase transition occurs at Pc = 1/2 [5]. Fisch 
[7], Schwartz |24] . and Aharony and Stephen [T] obtained equations to determine the transition 
temperature for square-lattice models with random interactions each of which takes Ji or with 
the same probability. 

The problem is simpler if the disorder or impurity exists only in one direction and the system 
maintains the translational invariance in the other direction. Let us consider an Ising model on the 
square lattice. Suppose that JL and denote coupling constants on a vertical and a horizontal 
bonds respectively. McCoy and Wu [5D] introduced a model where Jffs are constant and J^fs 
are uniform among the columns but row-to-row random. See Fig. [TJ Shankar and Murthy [25j 
considered a slightly different model, where jT’s are uniform among a row but row-to-row random 
while Jffs are uniform in the whole system (Fig. [T]). The criticality condition for both the models 
was derived in Refs. [liiin]. It is important that the partition function of the McCoy-Wu and the 
Shankar-Murthy models can be written as the largest eigenvalue, or the Lyapunov exponent in other 
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(a) McCoy-Wu (b) Shankar-Murthy 





Figure 1: Schematics for (a) the McCoy-Wu model and (b) the Shankar-Murthy model, (a) In 
the McCoy-Wu model, the interactions Jj^’s on the horizontal bonds are uniform and Jij's on the 
vertical bonds are random from row to row but uniform from column to column, (b) In the Shankar- 
Murthy model, J^j ’s are random from row to row but uniform among a row, while ’s are entirely 
uniform. 


words, of a product of random 2x2 matrices. As far as the transition temperature is concerned, 
however, it is sufficient to work with random diagonal matrices. 

In the present paper, we consider random Ising models with striped randomness on the hon¬ 
eycomb and square lattices. These models are extensions of the Shankar-Murthy model, and were 
discussed earlier by Hamm and Hoever |10) . The random couplings are not always ferromagnetic. 
As we shall see below, in these models one encounters computation of the Lyapunov exponent of a 
product of random non-diagonal matrices to determine the location of a phase transition. This is 
contrasted with the situation for similar but inequivalent models with ferromagnetic couplings dis¬ 
cussed in Refs. HU, for which the transition temperature is determined exactly by a single numerical 
equation. Unfortunately, there is no complete mathematical method to compute such a Lyapunov 
exponent. In the present paper, nevertheless, we successfully compute the exact Lyapunov exponent 
in the limiting case, where all the non-zero couplings have the same magnitude and the temperature 
is zero, and provide the exact location of the zero-temperature phase transition on the axis of the 
density of impurities. We furthermore provide precise phase diagrams on the basis of the numerical 
computation of Lyapunov exponents. 

The outline of the present paper is as follows. We first define the models in Sect. [2] We then 
introduce the transfer matrices and their Majorana-fermion representations in Sect. [S] We derive an 
equation which determines the transition temperature there. On the basis of the equation derived 
in Sect. |31 we visit known exact results for pure systems in Sect. SI We then give main results on 
random systems in Sect. 0 where the zero-temperature phase transition is first discussed and finite 
temperatures follow it. Section [B] is devoted to the conclusion. 

2 Models 

Let us begin with the honeycomb-lattice model. As shown in Fig. [U the honeycomb lattice is 
equivalent to the brick lattice. We write the Ising spin sitting on each lattice point of the brick 
lattice as Sij, where i and j label vertical and horizontal positions respectively. We then assign 
a spin-spin interaction J on the vertical bonds and Ji on the horizontal bonds, where i labels the 
vertical position of a bond. We assume that J^’s are random and independent for different i. The 
Hamiltonian is written as 

M N 

H = ( 1 ) 

i=i j=i 
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(a) honeycomb lattice 



(b) brick lattice 



(c) square lattice 



Figure 2: Lattices studied in the present paper, (a) The honeycomb lattice, (b) the brick lattice, and 
(c) the square lattice, (a) and (b) are equivalent. For the brick lattice (i.e., the honeycomb lattice), 
we assume that the interactions on the vertical bonds denoted by J are uniform, while those on 
the horizontal bonds denoted by Ji are uniform in a row but row-to-row random. As for the square 
lattice, we assume the same property of the interaction on the vertical bonds. The interactions on 
the horizontal bonds, however, are independently random in the adjoining two bonds and alternate 
in a row, and besides they are row-to-row random. 

MI'lNjl 

-EE iJ2fj.-lS2fi-l,2i^-lS2fj.-l,2v + J2fj.S2fj.,2ijS2fi,2v+l) , 

fl—1 ly—l 

where M and N denote the number of sites on the vertical and horizontal lines 
assume that M and N are even numbers. We consider two bimodal distributions 
follows: 

(±J model) P{Ji) =p6{Ji - Jq) -I- (1 -p)S{J^ + Jo), 

(diluted model) J’(Jj) = pd{Ji — Jq) -|- (1 — p) 6{Ji), 

where Jq > 0 and 0 < p < 1. Note that p stands for the density of ferromagnetic bonds in the 
horizontal bonds. We moreover restrict ourselves to the range of p such that [Jj] > 0, where the 
notation [• • •] means the average over the randomness of the bonds: 

[■■■]= I (3) 

i 

We next consider the square lattice. Same as the brick lattice, we assume that the interactions 
on the vertical bonds are uniform. Regarding the horizontal bonds, however, we assume that the 
interactions are random in the adjoining two bonds and form an alternating order, and besides they 
are independent in different rows. An instance of the configuration of interactions is depicted in 


respectively. We 
for random Ji as 

( 2 ) 
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Fig. He). The Hamiltonian of the present model is given by 


H 


M N 

i=l 0 = 1 
M N/2 

-EE {JilSi^ 2 v-lSi^ 2 v + Ji 2 Si^ 2 vSi^ 2 v+l) ■, 


i—\ v—1 


( 4 ) 


where Jn and Ja (z = 1, 2, • • •, M) are independently random and follow Eq. (H- 

For both lattices, we assume the periodic boundary condition in the horizontal direction {Si^N+i = 
Si^i) and in the vertical direction (^m+Ij = 

It should be noted that the ±J models on both lattices mentioned above are invariant under 
the change of the sign of all the horizontal couplings and that of the spins in every odd column. 
Therefore it turns out that the ± J models have the ferromagnetic-antiferromagnetic symmetry with 
respect to p = 1/2. 


3 Transfer matrices 


We first focus on the square-lattice model. A portion of the Boltzmann factor regarding the lowest 
two layers in Fig. is written as 


f N/2 N 

exp I {KiiSi^2v-lSl^2v + Ki2Si^2vSl^2v+l) + L SijS2,j 

i=i 


I !y=l 


N/2 N 

+ ''^^{K2iS2,2u-iS2,2u + K22S2,2uS2.2i^+i) + S2,jS3j j , 

i/=l j=l 


( 5 ) 


where Kij = Jij/T (z = 1, 2, • • •, M; j = 1,2) and L = J/T with the temperature T in the unit of 
fcs = 1. The transfer matrix between the lowest and the second lowest rows is given by 

TK^^,K^2,K2l,K■22{Sll^ • • •, S'iat; S'31, • • •, S^n) 

= E • • •)'5'lAr; «S'2i, • • •, <S'2Ar) 

S 21 -,---,S 2 N 

xAif2j_if22(S'2i, • • •, S2N] -S'31, • • •, S^n)^ (6) 


where 


AKii,Ki2{^iiN '' 5 ‘Siat; -S21, • • •, S2n) 

( N /2 N 

E (ArilS'i_2j/-lS'i^2jv + Ki 2 Si^ 2 uSl. 2 v+l) + L T.S 11 S 20 

v=l i=l 

Using these 2^ x 2^ transfer matrices, the partition function is written as 

Z = XT (TR-ii,_R-i2,_R'21,if22Fif3i,if32,if41.if42 ’ ’ ’ MM-1,2,^M1,KMi) ' 

The matrices in Eqs. 0 are now written as 

Aa-h^atijC-SiI) • • •, -Siat; S'21, • • •, S2n) = C{Sii, - ■ ■, Sin\Akm,Ki2 I'S2i, • • •, S2n) 


(n/2 


N 


A 


Kii,Ki2 


= exp E {Kii(T 2 u-i<^L + ^ 120 -Lo-L-ii) exp < 


\ 




( 7 ) 


( 8 ) 

( 9 ) 

( 10 ) 
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where cr“ {a = a;, y, z) is the a-component of the Pauli operator and • • •, <5'iiv)} is the basis of 
an iV-spin system satisfying ajlSu, • • •, Sin) = Sij\Sii, • • •, SiN) for j = 1, 2, • • •, TV. C and g are 
defined by 


C = ^isinh 25 ^ 


-W/2 


g = 2 logcothL. 


( 11 ) 


We introduce Majorana fermions and their Fourier transformations as follows; 




1 

71 


'2v-2 




cr 


z 

2u-l^ 




1 

71 



G 


y 

2y-l^ 


( 12 ) 


'2y-l \ 

n I ( 13 ) 

. k=l / 

= M E (cn(<z)e*^‘'+ct(g)e-«'^), (14) 

0<g<'7r 

where n = 1,2, 3,4 and ly runs from 1 to N/2. Although the possible values of the wave number 

q depend on the parity of 11^=1 they have no effect on the argument below and hence we do 

not care about them. Then Akii,Ki 2 in Eq. ([S]) is written as Akii,Ki 2 = ifi2)®^P(^77 

where 


■03 (i^) = 


72 


n Ml')= 




72 


H 


( 1 ) 

Kii,Ki2 


NI2 

2iE {KiiMi' - MAi') + Ki2MMi{i' + 1)) 
2i ^ \Kii {c2{q)cl{q) + c\{q)cAq)') 

0<q<7r 

+ Ki2 (e“*«C4(g)c|(g) + e*«c7g)ci(<3')) } 

0<g<7r 


(15) 


N/2 

= M E (^1 ( 9 ) 4(7 + 4 ( 9 )C 2 ( 9 ) + cAq)c\{q) + c\{q)ci{q)') 

0<<7<7r 

= E 4 '^( 7 - ( 16 ) 

0<g<7r 


Since each mode is decoupled from others, the transfer matrices are given by a product of those 
with a fixed mode. Hence one gets 


Z = 


It {^Kii,KM)M2i,K22{q) ■ ■ MKMi,KM2{q)'\ > 


0<(?<7r 


where 


= exp(h 7 i,Ki 77 )exp( 4 ^ 47 )- 

Taking the limit N ^ oo, the free energy per spin is written as 

where ^ 

/(*?) M ^"^^ 11,^12 (9)7if2i ,K 22 (.q)''' •^Kmi ,Km 2 ( 9 ) I' • 


(17) 


(18) 


(19) 

( 20 ) 
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In order to seek the critical temperature, we will investigate fo = limM->oo /(O)- We will see in fact 
that fo is non-analytic at a critical point Tc. 

In order to simplify AKii,Ki 2 {f^) further, we introduce new Majorana fermions as 


Cn{q) = {an{q) + iPn{q)) 


1 


Cn(<?)'^ = {0tn{q) - il3n{q)) ■ 

The transfer matrices are rewritten in terms of these new Majorana fermions as 

(a2a3 +/32/33) 

+2iKi2 [cosg (a 4 ai -I- PaPi) — sing ( 04 /?! -|- , 

(g) = 2ig (aia 2 + a 3 a 4 + Pih + / 33 / 34 ), 


( 21 ) 

( 22 ) 

(23) 


where we have omitted (g) on and Pn- It is clear in the above that and /3„ are completely 
decoupled and symmetric in the limit g —>■ 0. Hence the trace in Eq. (EHl) with g — ^ 0 is given by 
the square of the trace on the a subspace. 

Now, focusing on g —>■ 0, we transform Majorana operators into the Pauli matrices rf, 
and as 

ai = a2 = -^rf, a3 = “4 = (^4) 

Then we obtain 

2i{Kna2a3 + K^2aiai) = (ifn - ifi2H)TfT|, (25) 

2ig{aia2 +uoUa) = g{l + V)Tf, (26) 

where the unitary matrix V is defined as H = Tfrf. This matrix V commutes with 

(o') 

hg TO) and has eigenvalues ±1. Thus the transfer matrix is decomposed into two on the subspaces 
with V = ±1 respectively. Therefore /(O) in Eq. (1201) is written as 

, / M M \ 


/(O) = -2T — In ( tr A+ + tr A( 


(27) 


2=1 


2=1 


where Af and are 2x2 matrices on the subspaces with V = 1 and V = —1, respectively. We 
take the basis for the V = 1 sector as eigenstates of rfrl, i.e.. 


|+) = ^(|tt) + IU)), 


i-) = ^(in) + iit)). 


(28) 


Using this set of basis, Af is written in the matrix representation as 

/ 0 

* ^ 0 Q-iKii-Ki2) 

For the U = — 1 sector, we settle the basis as 

|+) = ^(|tt)-|U)), |-) = 

Then one obtains the matrix representation of A ~, 

n „-{Ki^+Ki2) 


cosh 2 g sinh 2 g \ 
sinh 2 g cosh 2 g J 

(29) 

^(in)-i;t)). 

(30) 


A- = 


(31) 


The second term in the brackets of Eq. (EH), *.e., the trace of a product of A^ is estimated as 
g 2 [Kij]M ]\/[^ since it is the trace of the diagonal matrices. Regarding the first term, we 

define the Lyapunov exponent of a product of random matrices by 


7 = lim — In tr TT Af. 

i=l 


M 


(32) 
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Using this, the first term is estimated as . Therefore fo is obtained finally as 

/o = -2Tmax{7,2[i^y]}. (33) 

The critical point is determined as a singular point of this quantity. 


3.1 Honeycomb lattice 

The model for the honeycomb lattice is obtained by making {Kii^Ki 2 ) = (0, Ki) for i = 2,4, • • •, M 
and {Kii,Ki 2 ) = {Ki, 0) for i = 1, 3, • • •, M — 1 in the above formulas on the square lattice with 
Ki = Ji/T. Therefore one gets the following result instead of Eq. (1551) . 

/o = -2Tmax{K,[/f,]}, (34) 


where 


= lim —Intr TT li?,.,, 

M^oo M 11 

/i=l 

n 


M/2 


= 


The relation between BJ and B- 


0 


aTKi 


cosh2g sinh 25 
sinh 25 cosh2g 


B- = UB+U, U = 


0 1 
1 0 


yields a simpler representation of Eq. (123 as follows. 


= lim ^IntrTTSj 

A/f —Vcvi A/7 


M 


M—yoo Ad 


where 


Bi = BTU = 


0 

0 


sinh2g cosh 25 
cosh 2g sinh 2g 


(35) 

(36) 

(37) 

(38) 

(39) 


4 Pure systems 

For the sake of confirmation, we visit the pure honeycomb-lattice model in the present section, 
before going to the random case. In the pure system, Ki is no longer random and, supposing 
Ki = K = Jq/T uniformly for any i, Eqs. (1^ . (l3^ and (l3^ reduce to 


fo = —2Tmax{K, K} 


and 

K= lim — IntrB^-^ = InA, 

M—¥00 A^I 

_ / 0 \ / sinh2g cosh2(; \ 

0 e~^ J cosh2(/ sinh2g J ’ 

where A is the largest eigenvalues of B. After a straightforward computation, one finds 


(40) 

(41) 

(42) 


A = cosh K sinh 2g + \J 1 + cosh^ K sinh^ 2g. 
Now, a simple algebra shows that the equation k = K reduces to 


tanhAT = 


1 

sinh 2L 


(43) 


(44) 
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We define the temperature as the solution of this equation. Since we have 


^ r -2 Jo for T < Tc 

■fo |-2TlnA forr>rc ’ 

T = Tc provides a singular point of /o and hence /. 

In the special case where L = K, namely J = Jq, Eq. (l4^ is further simplified into 


(45) 


cosh2i4r = 2. 


(46) 


This is nothing but the critical point of the Ising model on the isotropic honeycomb lattice m- 


5 Random systems 

The transition points are determined by the equation k = {Ki\ for the honeycomb lattice and 
7 = 2 [Kij] for the square lattice. In order to analyze these equations, we point out that the random 
matrices in Eq. (1^^ and Bi in Eq. (1551) have the same form as the transfer matrix of the one¬ 
dimensional random-field Ising model. Therefore the Lyapunov exponent k and 7 can be evaluated 
using techniques for the free energy of the one-dimensional systems. In this section, we obtain the 
exact value of the transition probability in the zero temperature limit and numerically calculate the 
phase boundary in the finite temperature, assuming J = Jq. 


5.1 Zero temperature 


In the zero temperature limit, the Lyapunov exponents 7 and k are regarded as ground-state energies 
of one-dimensional random Ising models. In general, it is a difficult task to obtain the exact ground- 
state energy of a random system. However, the exact Lyapunov exponent in the case of J = Jq can 
be obtained by the method proposed in Ref. [14]. First we focus on the honeycomb lattice with the 
±J distribution. 

We consider a one-dimensional Ising-spin system defined by the transfer matrices Bi in Eq. (1551) . 
Supposing that Zf denote the partition functions of this system with sites 1, 2, • • •, z under the fixed 
boundary condition S'i = ±1 respectively, they obey the following recursion relation 






= B.[ ), (47) 

"i+l 

where we assume the initial condition = Z^ = 1. The asymptotic behavior of Zf^ for z ^ 1 in 
the low-temperature limit K = J/T ^ 1 should be written in the form of 






(48) 


where z = and Xi and are random exponents to be investigated below. We here omitted 
prefactors, since they do not affect the propery of the exponent in the zero-temperature limit. The 
ground-state energy per spin of the one-dimensional system is expressed as 


K XM + max{aM,0} 

lim — = lim --- 

T-i-O K M->oo M 


1 “ 


(49) 


where we use the fact that the exponent om does not diverge with M as will be shown later. The 
last expression implies that the average increment of exponent Xi corresponds to the Lyapunov 
exponent. 

When Ji = -|-J, the random matrix Bi has the asymptotic form 


B, = 




(50) 



where we have used in Eq. (EH) 


sinh25~z cosh2(7^1. (51) 

Then the recursion relation of Eq. (|T7l) with Eq. (H51) is reduced to 

Xi+i = Xi — 1 + max{ai, —2} (52) 

Ci+i = max{ai, 2} — maxjai, —2}. (53) 

Similarly, when Ji = —J, we obtain 

Xi+i = + 1 + max{ai, —2} (54) 

Oi+i = max{ai, 2} — max{ai, —2} — 4. (55) 


Now, noting ai = 0, one can see that possible values of are restricted only to 0, ±2, and ±4. The 
transition matrix which transforms to a^+i is given by 



Ui 

= 4 

ai — 2 

ai =0 

ai — 2 

ai = 

-4 

<li+l = 4 

/ 

0 

0 

0 

P 

P 

\ 


ai+1 = 2 


0 

0 

P 

0 

0 



li+l = 0 


P 

p 

0 

P 

P 


(56) 

fli+i = —2 


0 

0 

P 

0 

0 



Oi+i = -4 

1 

p 

p 

0 

0 

0 

) 



where p = 1 — p. Since the Markov chain generated by this transition matrix is clearly irreducible 
and aperiodic except for non-random cases, it is ergotic and has the unique stationary distribution 
corresponding to the maximum eigenvalue of one. The distribution of ai converges in the limit 
i > 00 to this stationary distribution which is explicitly given by 


^a(+4) 

Pa{0) = 


Pa{-2) 


PP 

1 -l-p’ 


Pa{ + 2) 


l-pp 


(1 + P )(1 + P )’ 

^ P(1 - PP) 

(1 + P )(1 + P )’ 


P{^-PP) 

{l+p){l+p)’ 


Pa{-i) 


PP 

1 -l-p' 


The increment of the exponent xt depends on Ji and max{ai,—2} as shown Eqs. (15^ and (15^ . 
Therefore, its average increment is computed by 


{-p + p) + {4Pa(+4) + 2P,(+2) - 2Pa{-2) - 2P„(-4)} , (57) 


which yields the exact Lyapunov exponent in the zero-temperature limit 

lim^= 


T^O K {l+p){2-p)' 


(58) 


Finally, noting \Ki] = {2p — \)K for the ± J model, we obtain the exact critical probability Pc 
in the zero-temperature limit as 


Pc = 1 — 2 sin — = 0.65270364 • • •. 
18 


A similar calculation for the square lattice yields 


and the critical probability 


lim^ = 2p(l-p), 


\/5- 1 

Pc = — = 0.61803399 • • •, 


(59) 

(60) 

(61) 
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Figure 3: p dependence of the Lyapunov exponent for the ± J model on the honeycomb lattice. The 
dashed line indicates [Ki\/K. The critical point corresponds to an intersection point between the 
curve of the Lyapunov exponent and the dashed line. 


which is equal to the inverse of the golden ratio. 

In the diluted models, the Lyapunov exponents are obtained as 


ito “ = hijih. 

T-).o K 3 - p 


ii„ 1 = 

T-s-o K 3 


(62) 


The critical probability in the zero-temperature limit is Pc = 0 both for the honeycomb and square 
lattices. This result as well as the known fact that the pure system (p = 1) has the ferromagnetic 
ground state lead to the conclusion that the zero-temperature phase of the diluted model for p > 0 is 
ferromagnetic. This conclusion is naturally understood since this model is always percolated except 

p = 0. 


5.2 Finite temperatures 

So far some methods have been developed to numerically compute the Lyapunov exponent of a 
product of random matrices. Mainieri proposed the cycle expansion method whose convergence is 
exponentially fast in the cycle length m- Bai improved this method using the evolution operator 
approach [5] . We emply this improved cycle expansion (ICE) method and the Monte Carlo method. 
The ICE method is summarized briefly in Appendix. Although we still consider J = Jq for com¬ 
parison with the zero-temperature results, numerical methods used in this section can be applied 
to J ^ Jq. 

For the honeycomb-lattice model with T/ J > 0.7, we used the ICE with the cycle length n = 20, 
which produces more than 15 reliable digits. Eor lower temperatures, however, the ICE needs a 
still longer cycle length to attain the convergence. This cost was complemented with the Monte 
Carlo method for T/J< 0.7, where Lyapunov exponents were obtained by averaging over 1000 
configurations of a product of 10^ matrices with the standard error of n/K less than 6.0 x 10“^. 

As for the square lattice, since the random matrix has three possible choices, the ICE needs 
more computational cost as well as longer cycle length for convergence than the honeycomb lattice. 
Therefore we used the ICE with the cycle length n = 15 for T j J > 1.0 and the Monte Carlo method 
for T j J < 1.0. Results by ICE have the same accuracy as the honeycomb lattice, while the standard 
errors of q/AT by the Monte Carlo method are less than 9.0 x 10“®. 
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Figure 4: p dependence of the Lyapunov exponent for the diluted model on the honeycomb lattice. 


The p dependence of the Lyapunov exponent for the ±J model on the honeycomb lattice is 
shown in Fig. [31 The solid curve indicates the exact result (j58|) in the zero-temperature limit. The 
critical points are determined by the intersection point with the dashed line {\Ki]/K = 2p— 1). The 
zero-temperature critical point is given exactly by Eq. (1331) . In Fig. |T]we show the results for the 
diluted model. It is clear that the critical probability in the zero-temperature limit is equal to zero 
as mentioned before. 

The critical line in ihe p — T plane were obtained by solving the equation k = [Ki] or 7 = 2[Kij\. 
For the temperature region where the ICE is valid, we estimated the critical temperature using 
the bisection method for given p’s. On the other hand, for lower temperatures, we obtained the 
critical probability for given T’s using the Monte Carlo method and the bisection method. The 
results are shown in Eig. [3] Although one cannot say anything about the property of the phases 
away from the critical line, it has been shown exactly that the critical temperature at p = 1 is 
the transition point between the high-temperature paramagnetic phase and the low-temperature 
ferromagnetic phase. This fact as well as the fact that the critical line obtained here is continuously 
connected with the p = 1 transition point imply that the critical line is in fact the transition line 
separating the paramagnetic and ferromagnetic phases. It is worth noting that the ± J model has 
the ferromagntic-antiferromagnetic symmetry with respect to p = 1/2, as mentioned in Sect. [31 
Therefore, we deduce that there is a striped antiferromagnetic phase in the small p region, though 
it is not shown in Fig. 0 In addition, the study on the Shanker-Murthy model suggests that the 
Griffiths phase may exist between the ferromagnetic and antiferromagnetic phases and below the 
critical temperature of the pure system [51 [33] . 

We comment that Hoever nni has numerically studied the critical line of the ± J model on the 
square lattice. Our result on the same model is in perfect agreement with it. 

It should be pointed that the critical lines of the ± J model on both lattices are not vertical near 
zero temperature on the p—T plane. As shown in Fig. [31 the Lyapunov exponent k devided by K 
at low temperature seems to have a form 

k/AT ~ eo(p) -I- cT, (63) 

where eo(p) denotes the zero-temperature limit of n/K which is exactly given by Eq. (I58|) . The 
coefficient c is estimated as c = 0.1509 by fitting the data at p = pc for 0 < T/J < 0.2. This 
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Figure 5: Critical lines in the p—T plane of the ± J models and the diluted models on the honeycomb 
and squara lattices. 

observation implies that the critical temperature for T/J <C 1 is behaves as 

Tc{p) ^ {(2p - 1) - eo(p)} ~ 15.70 x (p - p^), (64) 

which is consistent with our numerical result (Fig. [^. On the square lattice, we obtain Tc(p) ~ 
14.64 X (p- Pc)- 

These behaviors of the critical temperature are quite different from the Shankar-Murthy model, 
where the constraint Kn = Ki 2 leads to the exact Lyapunov exponent 7 = 2p regardless of p. Then 
an essential singularity of p ^ at T = 0 causes the vertical growth of the critical temperature 

at Pc = 1/2. 

In contrast to the ± J models, the diluted model has the vertical phase boundary in the same way 
as the Shankar-Murthy model. The Lyapunov exponent at p = 0 is exactly calculated as k = 2p 
for the honeycomb lattice, which is essentially singular with respect to T. The low-temperature 
behavior of the Lyapunov exponent is shown in Fig. [T] 

6 Conclusion 

In this paper, we discussed the two-dimensional Ising model with striped randomness on the honey¬ 
comb lattice and on the square lattice. We simplified the transfer matrices by using the Majorana 
fermion operators and obtained the 2x2 matrix representation in the long-wavelength limit. The 
Lyapunov exponents can be calculated by regarding these decomposed matrices as the transfer 
matrix of the one-dimensional random-field Ising model. We obtained its exact solution in the zero- 
temperature limit and determined the exact value for the critical probability Pc of ferromagnetic 
bonds. The critical line on the probability-temperature (p — T) plane was also calculated with highly 
accurate numerical methods. 

We focus only on the zero wave-number limit to analyze the critical point. However the wave- 
number dependence of the free energy is necessary to analyze critical phenomena and to determine 
the universality class, which remains as one of the future issues. 
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Figure 6: Low-temperature behavior of the Lyapunov exponent for the ± J model on the honeycomb 
lattice. The solid line is obtained by fitting the data at p = pc for 0 < T/J < 0.2. The intersection 
with the horizontal dashed line indicates the transition temperature. 



Figure 7: Low-temperature behavior of the Lyapunov exponent for the diluted model on the hon¬ 
eycomb lattice. 


13 










A Cycle expansion of the Lyapunov exponent 

In this appendix, we briefly summarize the improved cycle expansion method as a numerical method 
to compute the Lyapunov exponent. 

The Lyapunov exponent of a product of random matrices is defined as 

7= lim — (In IIM1M2 • • • M„||) (65) 

n—>-oo 77, 

where Mi is a random matrix and the angular brackets denote the ensemble average. It is no¬ 
table that the Lyapunov exponent is independent of a choice of the matrix norm || • || as long as 
equivalent norm used. The maximum eigenvalue /io(M) is the most convenient norm to derive the 
cycle expansion because it is invariant under cyclic rotation /io(MiMj) = fj,o{MjMi) and satisfies 

Mo(M2)=/ro(M)2. 

Assume that the random matrix Mi has m possible choices to be Mi = with probability Ps 
(s = 1, 2, • • •, m). With a string of n choices S = S1S2 • • • Snj the product of n random matrices 
As^As^ ■ ■ ■ Ag^ is denoted by As and its probability is by P{S) = W^^iPsi- With these notations, 
the Lyapunov exponent is expressed as the large-n limit of the following quantity 

7n = - V] P(5')ln/xo(As), (66) 

n ^' 

|S|=n 


where |S'| is the length of a string S. 

Owing to the properties of /ro(A), the sum in Eq. (IHHl) is decomposed into the sum over the set Cp 
of all possible primitive cycles, where a cycle is said primitive if it is not a repeat of a smaller length 
cycle. Two cycles are equivalent if they differ only by a cyclic rotation. For example, S1S2S1S2 ^ Cp 
because it is a repeat of S1S2, and if S1S2S3 G Cp, S2S3S1 ^ Cp. As a result, the cycle expansion of 
the Lyapunov exponent is obtained as 

n 

7™=EE <5,|S|,„P(5)nn/xo(As). (67) 

SeCp r=l 


This formula was also derived from the expansion of the thermodynamic zeta function and its 
exponential convergence with the cycle length was observed m- 

Bai proposed an accelerated algorithm for the cycle expansion based on the evolution operator 
approach [2] . Bai numerically showed super-exponential convergence of the weighted average of the 
cycle expansion. 


In = 



( 68 ) 


in contrast to exponential convergence of the cycle expansion method. The averaging weight Wk is 
determined to eliminate all exponential converging terms by using the evolution operator approach. 
In the case of 2 x 2 matrices with positive elements, Wk is given by 


Wk = - 


1 

k 


k 

t=l 


(69) 


I c-i 

= 2^ — 


l-o(As)'- ’ 

SeCpr=l 

where g{A) denotes the ratio of the second eigenvalue of A to the first one po{A). 


(70) 
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